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ABSTRACT 

A  finite  element  method  for  structural  optimiizaticr.  of 
a  prism.atic  beam,  of  a  hom.ogeneous  ,  isotropic  material  is 
developed.   The  beam  has  a  rectangular  cross-section  c: 
constant  fixed  height,  a  fixed  length,  and  a  fixed  vol--.,r-:e. 
Structural  optim.um  is  defined  as  that  shape  which  allcv-  a 
maximum  load  within  the  elastic  range. 

A  computer  program  is  developed  to  solve  the  resulTing 
system  of  equations  and  various  example  problems  are  solved 
Comparison  is  m.ade  with  exact  optim.um  beam  designs  v/here 
possible. 

The  finite  element  model  is  able  to  solve  problem.s  with 
any  boundary  conditions  and  types  of  loading  that  are  : in- 
sistent with  the  number  of  elements  selected. 
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I.   INTRODUCTION 

Optimization  has  a  history  as  old  as  the  universe  itself. 
The  planets  have  taken  their  optimum  position  in  the  solar 
system,  the  solar  system  in  the  galaxy,  etc.   The  whole 
history  of  evolution  is  one  of  gradual  optimization.   Man 
has  used  optimization  in  one  form  or  another  since  the 
earliest  of  times.   Individual  man  was  weak  with  respect  to 
his  hostile  environment.   He  learned  to  become  a  social 
animal  so  that  as  a  group  he  was  strong  -  optimization. 
In  wars  between  groups  spears  defeated  rocks,  arrows  defeated 
spears,  and  so  on  to  the  present  day.   All  were  essentially 
optimization . 

The  early  Egyptians  formalized  some  optimization  tech- 
niques in  their  development  of  geometry  and  trigonom.etry . 
The  developm.ent  of  the  Calculus  enabled  us  to  find  optimum 
function  values.  The  Calculus  of  Variations  made  possible 
the  optimization  of  functionals,  the  optim.um  function  from 
a  function  space. 

In  structural  optimization,  however,  oftentimes  the 
Calculus  of  Variations  produces  nonlinear  differential 
equations  that  cannot  be  solved  in  closed  form  and  the 
numerical  methods  of  solution  are  difficult.   If  these 
differential  equations  could  be  transform.ed  into  algebraic 
equations,  a  solution  might  be  more  readily  obtained. 


The  finite  element  rr^ethcd  has  been  applied  to  the 
theories  of  elasticity  [1],  [2]  and  mechanical  vibrations 
[3] J  as  well  as  to  other  fields,  to  resolve  such  problems. 
It  has  been  used  in  the  optimization  of  truss  networks  but 
as  yet  not  in  actual  beam,  optimization.   Since  the  m.ethod 
is  essentially  variational  in  nature,  it  should  have  an 
application  to  structural  optimization  theory.   This  sug- 
gests the  desirability  for  research  in  this  field. 

Since  this  is  a  first  attempt  at  a  formulation  of  this 
type,  we  v/ill  restrict  ourselves  to  a  specific  type  of 
problem;  the  structural  optimization  of  a  prismatic  beam 
of  one  material  which  is  homogeneous  and  isotropic.   This 
beam  has  a  rectangular  cross-section  of  uniform  fixed  height, 
and  has  a  given  volume  and  length.   The  structural  optimum 
in  this  case  is  the  shape  required  to  maximize  the  load 
within  the  elastic  range. 

A  variational  formulation  is  first  presented  to  establish 
the  optimization  technique  for  isoperim.etric  problems.   The 
finite  element  method  is  applied  to  the  problem  using  a 
parallel  approach.   A  model  is  developed  and  then  checked 
for  validity  by  solving  various  applicable  problems. 
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II.   VARIATIONAL  FORMULATION  OF  BEAM  OPTIMIZATION  PROBLEMS 

A.   DEFINITION  OF  THE  PROBLEM 

Salinas  [^]  shows  that  the  problem  of  structural  opti- 
mization of  a  beam,  under  a  fixed  volum.e  constraint,  m.ay  be 
solved  by  using  calculus  of  variations  techniques.   He  shows 
that  operation  with  the  techniques  for  isoperimetric  prob- 
lems [5]  on  the  augmented  potential  energy  functional, 
T*  =  T-XS,  produces  equations  which  describe  the  struc- 
turally optimum  beam.   T  is  the  potential  energy  of  the  beam 
and  S  is  the  isoperim.etric  constraint.   The  quantity  X  is 
a  constant  Lagrange  multiplier  which  is  actually  a  measure 
of  the  strength  of  the  beam. 

We  shall  consider  the  problem  of  structural  optim.ization 
of  a  sim.ply  supported  beam  of  length,  L,  and  volume,  V  . 
The  applied  load  is  uniform,  P   (force/unit  length).   The 

cross-section  of  the  beam  is  such  that  the  m.oment  of 

n 
inertia  may  be  given  as  I  =  CA   where  C  and  n  are  constants. 

Gravity  forces  and  shear  effects  are  assumed  negligible  and 

the  usual  beam  assumptions  [6]  are  used. 

i)  sm.all  deflection  theory 

ii)  uniaxial  stress  state 

Figure  1  graphically  depicts  the  beam  and  defines  the  x 

and  y  axes.   VJe  seek  the  shape  of  the  beam,  to  give  a  maximum 

load  intensity  subject  to  the  constant  volume  constraint. 
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A(x)  (not  necessarily  rectanpular' 


y  ,v 


Figure  1.   Simply  Supported  Beam, 


B.   THE  VARIATIONAL  FORMULATION 

We  first  form  the  potential  energy  functional  from 
strength  of  materials  considerations.   The  displacement  of 
the  beam  Is  denoted  by  v  and  the  primes  Indicate  differen- 
tiation with  respect  to  x. 


L 


T  =  / 
0 


[¥ 


(v")^  -  P  V 
o 


dx 


The  fixed  volume  constraint  m.ay  be  written  as  /  A(x)dx=V 

0         ^ 
thus  the  constraint  augmentation  becomes: 


As  =  Xf      A(x)  dx 
0 


n 


Noting  that  I(x)  may  be  denoted  as  CA(x)  ,  the  augmented 
potential  energy  functional  Is: 
L  r^^„  ,__Nn 


T«  =  / 

0  L 


^^^^^^   (v")2  -  p  V  -  AA(x) 
2  o 


dx 


where  v  and  v"  are  v(x)  and  v"(x)  the  displacement  and 
curvature  of  the  beam  respectively. 
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In  order  to  extremize  the  above  functional,  we  first 
take  the  variation  [7]  of  the  functional  v/ith  respect  to 
displacement,  v,  and  set  this  equal  to  zero.   The  result 
is  the  differential  equation  of  equilibrium.   This  variation 
yields : 

(EAC^v")   -  P   =0 
o 

which  may  be  written: 

(a'^v")"  =  P  /EC  . 

Integrating  this  equation  twice  with  respect  to  x  gives: 

/  2  \ 

A^v"  =  P  /EX   —  +  K,  X  +  K„ 
o     I  2     1      2  1 

The  constants  of  integration,  K,  and  Kp ,  are  obtained 
from  the  boundary  conditions.   For  a  simply  supported  beam, 
the  end  moments  are  zero,  or: 

(a)  A^(0)v"(0)  =  0 

(b)  a"(L)v"(L)  =  0 

Boundary  conditions  (a)  and  (b)  yield  K^  =  0,  K,  =  -L/2 
respectively . 

Thus  we  obtain  the  differential  equation  of  equilibrium 
for  the  simply  supported  beam. 

A^v"  =  2E§  ^>^^-Lx)  ^"-1^ 

Taking  the  variation  of  the  functional  with  respect  to 
area,  A(x),  and  setting  this  equal  to  zero,  gives  the  dif- 
ferential equation  which  we  refer  to  as  the  optimality 
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condition.   For  our  example  the  optirality  condition  is 
HEC  (;,_)n-l(^„^2  _  ^  ^  0 

Solving  for  v"  yields: 

'n-1 


v"  =   /  2A   A  '  "  '  (II-2) 

v/nEC 


n-1 


Equation  (II-2a)  shows  that,  for  an  optirrium  beam,  the 
curvature  is  constant  only  for  n  =  1.   Otherwise  the  product 
of  a  power  of  the  cross-sectional  area  times  the  curvature 
is  constant. 

Multiplying  both  sj.des  of  equation  (II~2)  by  A   and  sub- 
stituting into  equation  (ll-l),  the  development  continues  as: 

n+1 
From  equation  (II-l): 

nil  , p 

.  2   fTT^      o  .  2  T  s 

Hence, 

nP 
A n+1      o   /  2  T  N 2 

A    =  BecX  ^^  -L^^ 

2     2         2  2 
Since  (x  -Lx)   =  (Lx-x  )   it  is  convenient  in  further 

development  to  use  the  latter  expression.   Substituting  this 

into  the  previous  equation  yeilds: 


1^4 


2 
„n+l      o  ,,         2v^ 

Taking  the  (n+D^h  ^oot  of  both  sides 

_1 
A  = 


■"nP^        2  2" 


n+1 


(11-^) 

Introducing  the  constraint  equation  and  solving  for  \   v;e 
obtain : 

1 


L  L 

/  A(x)dx  =  / 
0  0 


'"nP^        2 
BIcT  "^Lx-x-) 


n+1 


dx  =  V 


o 


But  since  ^  is  a  constant  we  obtain, 

1_ 

\    ^ 
0 


A  = 


^2\n+l 
nP,  \ 


2 


n+1 


Hec  I 


V 


o  / 


2  ^+1 
/  (Lx-x  )    dx 


(II-5) 


Rearranging  (II-5)  and  solving  for  P   v/e  obtain  the 
equation  defining  the  maximum  load  intensity 


P   = 
o 


n+1 

^8ECr(v  )  2 
n     o 


2 


f     /T-    2vn+l  , 
/  (Lx-x  )     dx 

0 


n+1 


(II-6) 


Substituting  (11-5)  into  (11-^)  we  obtain  for  the  opti- 
mum area: 

2 


A  = 


,r  /T    2sn+l 
V  (Lx-x  ) 
o 

2 


(II-7) 


r    fy         2  xn+1  , 
/  (Lx-x  )    dx 

0 
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It  should  be  noted  here  that  the  integral  is  not  defined 
in  closed  form  for  all  n  and  therefore  must  be  numerically 
approxim.ated  for  som.e  problems. 

Summarizing  the  equations  that  define  the  solution  to  a 
given  problem  we  have: 

2 


A  = 


V  (Lx-x  ) 
o 

2 


(II-8) 


/  (Lx-x  )    dx 


P   = 
o 


n+1 
^BECX  .^j    s    2 
n      o 


2 


r    /T    2xn+l-, 
/  (Lx-x  )    dx 

0 


n+1 


(11-9) 


The  value  of  X  is  determined  as  follows  .   Fromi  equation 
(ll-2a) : 


2A 


nECA 


n-1 


From  strength   of  m.aterials    theory 


M  =    -   E(Iv") 

max  max 


If  we  define  the  yield  stress  of  the  beam  material,  S  ,  as 
the  maximum  allowable  stress,  we  obtain: 


M     R 
'  max  m.ax 


where  R    is  the  vertical  distance  to  the  furthest  fiber 
max 

from  the  neutral  axis .   Thus : 
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E(Tv")    R 

Q   _  'rax  r.ax    ^  „    „ 

S   =  =  Ev     R 

y         I  max  r.ax 


Substituting  for  v" : 


S   =  ER     /   2X 


y      "^^V  ECA^-^ 


max 


^2  ^  2AE  I   R^ 


y     n   l^^n-l 


max 


nS   /„.n-l\ 

\   ^   /max 

It  may  be  shown  by  exam.ple  that  the  quantity  ^r-   is 

r2 

independent  of  x  and  thus  that  X    is  indeed  a  constant  for  a 
given  type  of  cross-section. 

The  first  example  is  a  constant  height  rectangular  cross- 
section  of  variable  width.   The  height  and  width  are  h,  and 
b(x)  respectively. 

I(x)  =  ^iflll^=  b(x)h.^ 

A(x)  =  b(x)h 

,2 
I(x)  =  ^  A(x) 

,2 
.*.   C  =  ~  ,   n  =  1,   R  =  h/2 


Ca"  -'-  ^   _C_  ^   \? /\2   ^    1 
R^     r2    \?/\\  3 
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Thus  for  a  rectangular  cross-section  as  above,  equation 
(11-10 )  becomes : 
,2 
6E 


S' 


(11-11) 


The  second  example  is  a  circular  cross-section  of 
radius  r (x) . 


Kx)  =  I£W_=  (.r{x)2)'  ■  (i) 


A(x)  =  Trr(x) 


I(x)  =  -p-  A(x)^ 


C  =  ^  ,   n  =  2,   R  =  r 


Ca"  "^    CA    /I 


r2     r2    ^^- 


(irr  )  • 


r2, 


Thus  for  a  circular  cross-section,  equation  (11-10 
becomes : 

S  2 
X    =   tI-  (11-12) 


4E 


The  third  example  is  a  rectangular  cross-section  of  con- 
stant width,  b,  and  variable  height,  h(x). 


!(,)  =  !M|)i=  b3h(x)3 


A(x)  =  b  •  h(x) 

I(x)  =  — ^  A(x)3 
12b 


1 


12b' 
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/.   C  =  — ;t  ,   n  =  3,   R  =  h/2 

Ca"~^    CA''     1     /k2v,2v    I  ^  \        1 

2~   "  ~?~  "  2  ^b  h  )  •   —  =  - 

R     R      12b  Ih  /    ^ 

Thus  for  a  rectangular  cross-section  of  constant  width 
and  variable  height,  equation  (11-10)  becomes: 

S  2 
A  =  5^  (11-13) 

Note  that  equations  (11-11),  (11-12),  (11-13)  depend 

only  on  cross  section  type  and  not  on  the  manner  of  loading 

L 
(independent  of  the  term  /   P(x)v(x)dx).   They  are  also 

0 

independent  of  the  boundary  conditions  since  they  were 
derived  from  equations  in  which  the  boundary  conditions  had 
not  been  envoked.   This  observation  will  be  used  in  Section 
111. 

Equations  (11-8),  (11-9),  and  (11-10),  when  applied  to 
a  specific  problem,  yield  the  structurally  optimum  beam 
shape  for  maximum  load  intensity  for  beams  which  conform  to 
the  assumptions  made  at  the  beginning  of  this  section.   These 
beams  must  be  simply  supported  and  under  a  uniform  load 
distribution. 

Appendix  A  gives  three  examples  of  structural  optimiza- 
tion problems  solved  with  equations  (11-8,  9,  10). 
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Ill .   A  FINITE  ELEMENT  MODEL  FOR  BEAM  OPTIMIZATinN 

A.  SPECIFICATIOIIS  AND  ASSUMPTIONS 

This  development  parallels  that  of  the  continuous  vari- 
ational fornulation  of  the  preceding  section.   V^e  consider 
prismatic  beams  of  one  m.aterial  which  is  homogeneous  and 
isotropic.   The  modulus  of  elasticity  and  yield  stress  of 
the  material  are  E  and  S   respectively.   The  cross-sectional 
area  of  the  beam  is  rectangular  with  constant  height,  h,  and 
variable  width,  b(x).   The  length  and  volum.e  of  the  beam  are 
denoted  by  L  and  V   respectively. 

The  plane  of  loading  is  in  the  centroidal  plane  parallel 
to  the  y  axis.   Gravity  forces-  and  shear  effects  are  assumed 
to  be  negligible.   The  usual  beam  assumptions  [6]  are  taken: 
i)  sm.all  deflection  theory 
ii)  uniaxial  stress  state 

B.  STATEMENT  OF  THE  PROBLEM 

As  in  Section  II,  we  seek  to  optimize  the  shape  of  the 
beam  (cross-sectional  area  as  a  function  of  position  along 
the  beam  axis)  to  give  a  maximum  load  intensity  subject  to 
the  constraint  that  the  beam  volume  is  fixed.   The  finite 
element  method  results  in  an  approximate  solution  to  this 
problem. 
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C.  DEFINITION  OF  THE  FINITE  ELEMENT 

The  beam  of  Figure  1,  v/lthout  supports,  is  discretized 
into  finite  elements  of  equal  length  as  shown  in  Figure  2. 
The  elem.ents  are  identified  by  num,bering  them  sequentially 
from  the  origin.   The  x  (centroidal)  axis  and  y  axis  (posi- 
tive downv/ard)  define  the  plane  of  loading  and  bending. 

V/e  consider  a  typical  element,  say  the  i   ,  as  shov/n  in 
Figure  3-   The  displacements  and  slopes  at  the  left  and 
right  ends  of  the  element  are  denoted  by  v,  ,0_,  Vp ,  9 
respectively.   The  cross-sectional  area  is  A..   The  local 
coordinate  system  x.  and  y.  are  defined  for  each  element. 

D.  APPLICATION  OF  THE  FINITE  ELEMENT  METHOD 

In  accordance  with  the  finite  element  method  [1] ,  shape 
functions  are  assumed  for  the  displacement  of  each  elem.ent. 
Since  a  cubic  polynomial  function  for  v.,  the  displacement 
over  the  i    element,  yields  continuity  of  displacements 
and  slopes  at  each  elem.ent  boundary,  it  is  taken  as  the 
assum.ed  displacement  field.   Thus  the  displacement  over  the 
i    element  m.ay  be  written: 

v(x. )  =  <N(x.  )>  {v} . 
1         1       1 

1  X  H        H  xl 

where  <N(x.)>  is  the  cubic  shape  function  vector  as  developed 
in  Appendix  B-1  and  {v}.  is  the  colum.n  vector  (v.,  9.,  v^,Q^) 


1 


The  elemental  cross-sectional  area.  A.,  is  assum.ed  to  be 

'   1 ' 

constant  over  the  element  to  obtain  a  tractable  mathem.atical 
m.odel.   This  point  is  further  discussed  later  in  the 
developm.ent . 
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y,v 


*—  X 


Figure  2.   A  Finite  Element  Beam. 


«*  X 


Figure  3.   The  I    Finite  Element. 
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since  the  augmented  potential  energy  functional, 
T*  =  U-W-AS,  must  be  extremlzed,  finite  element  formulations 
must  be  obtained  for  the: 

1)  strain  energy  of  bending  (U) 
11)  potential  energy  of  external  force  (W) 
ill)  Isoperim.etrlc  constraint  (  S) 

1.   Strain  Energy  of  Bending 

The  strain  energy  of  bending  of  a  beam  may  be  given 
as  : 

U  =  /   ^  (v")^  dx 
0   "^ 

The  finite  element  assumption  that  each  element  is 

th 
under  uniaxial  stress  enables  the  strain  energy  in  the  1 

element,  u. ,  to  be  expressed  as: 

A  EI.      ^ 

u.  =  /   — -    (vV)'^  dx. 
1    Q    2     1 

where  vV  denotes  the  second  derivative  of  v(x.)  with  resDect 
1  1  ^ 


to  X. 

1 


For  a  rectangular  section  of  constant  height,  h,  and 


width,  b. ,  the  moment  of  inertia  of  any  section  of  the 
element  may  be  written  as: 

I.  =  CA.  C  =  h^/12 

1      1 

The  strain  energy  in  the  1    element  is  obtained  as 


follows : 


v(x. )  =  <N(x. )>  {v}. 
1         1       1 
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v"(x..  )    =    <N"(x.  )>    {v}. 

(v")^    =    {v}T    <M">'y    <M">.     {v}. 

u.    =    /      ^  A.     {v}:    <N">T    <N">.     {v}.    dx. 
ip^2i  1  1  111 

1x4      Hxl        1x4        4x1 

Since    the   product    of   the    second   derivative    of   the 

T 
shape    function    vector,  <N'.'  N'.'>,  is    the    only    function    of 

x.    in   the    expression,    the    remaining   terms    may    be    brought 

outside    the    integral,    and   the    integral   may   be    evaluated. 

Thus  . 

u.    =    i^  A.     {v}:    /      <M">!    <N">.    dx.     {v}. 
1  2i1q  1  1         11 

The  integral  is  evaluated  in  Aj)pendix  B-2  and  called  [k*^], 
a  4x4  matrix.   Therefore: 

u.  =  ^  A.  {v}^  [k«]  {v}. 
1x4  4x4   4x1 

Here  it  may  be  noted  that,  if  the  A.  has  a  shape  function 

vector  associated  vilth   it  other  than  a  constant,  the  vector 

T 
{v}.  could  not  be  removed  from  the  integrand.   Since  the 

T 
values  for  the  vector  {v} .  are  unknown,  this  vector  would 

appear  in  the  integrand  between  tv;o  shape  function  vectors  , 

and  therefore  closed  form  integration  is  not  possible. 

When  the  term  EGA.  is  included  in  the  integrand,  the 

1 

result  of  integration  is  comm.only  known  as  the  elem.ent 
stiffness  matrix.   In  the  present  formulation  [k^]  is  called 
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the  "modified  element  stiffness  matrix"  since  the  kv .  are 
m.easures  of  elem^ent  stiffness. 

Since  strain  energy  Is  a  scalar,  the  total  strain 
energy  In  the  beam  Is  the  sum  of  the  Individual  element 
strain  energies.   Hence  v;e  have: 

U  =   Z   u.  =   E   ^  A.  {v}!  [k«]  {v}. 
1=1      1=1 

where  n  Is  the  number  of  elements. 

2 .   Potential  Energy  of  External  Forces 

The  potential  energy  of  external  forces  for  the  1 
element  may  be  written  as : 

w.  =  P   <D>.  {v}. 

10      11 

1x4   4x1 

The  vector  <D.>  Is  called  the  elem.ent  consistent  load  dls- 
1 

trlbutlon  vector  and  Is  described  In  detail  In  Appendix  B-3 

Since  the  potential  energy  Is  also  a  scalar,  the 
total  Is  the  sum  of  the  contributions  of  the  Individual 
elements.   Therefore: 


n 

W  =  P    Z   <D>.  {v} . 
o  .  .      11 
1  =  1 

3 .   The  Augmented  Potential  Energy  Functional 

The  constant  volum.e  constraint  In  finite  element 

n 
form  Is   S   A.A  =  V  .   In  order  to  maximize  the  load  Inten- 
1  =  1   1     o 

slty  P   we  form  the  augmented  potential  energy  functional, 

n  rn  n  n 

T*  =   Z   —  A.  {v}:  rk*l  {v}.  -   Z   P   <D>.  {v}.  -XI      A.A 
1  =  1  2   1     1  ^   ^     1    .^^   o     1     1     .^;^   1 
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E.   THE  FINITE  ELEP^.ENT  VARIATIONAL  PROBLEM 

The  extreir.izatlon  of  the  potential  energy  functional 

must  be  done  globally,  i.e.,  over  the  entire  beam.   A  tranS' 

formation  is  therefore  made  from  the  local  coordinates  {v}. 

to  a  global  coordinate  system,  {v.}.   The  transformation  is 

made  by  numbering  the  components  of  {v.}  sequentially  from 

the  origin,  v/here  the  first  component  is  the  displacem.ent 

at  the  origin,  the  second  Is  the  slope  at  the  origin,  etc. 

Thus: 

1 


^1 

■" 

^1 

^2 

= 

^l 

^3 

= 

1 
"2   = 

2 

^D 

= 

^l- 

< 

Therefore  the  global  representation  of  the  local  "dis- 
placement" vector  {v} .  becomes: 


{v} 


^^2i-l^ 


v 


V 


2i 


2i  +  l 


V^21+2^ 


(v.) 


In  global  coordinates  the  potential  energy  functional 

becomes : 

FT   "         T  n        _        n 

T*  =  ^   E   A.  {v.}^  [k«]  {v.}  -  P    E  <D>.  (v.)  -  XA  E  A. 
11             1       O.-i      11        ^_-,i 


2 


i  =  l 


i  =  l 


i  =  l 
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The  algebraic  equations  describing  the  structurally 
optimum  beam,  shape  and  maxim.um  perm.issible  load  intensity 
arise  from  taking  the  partial  derivatives  of  the  global 
augm.ented  potential  energy  functional  ivith  respect  to  the 
state  variables  of  individual  slopes  and  displacem_ents  and 
also  with  respect  to  the  control  variables  of  the  elem.ent 
areas  and  setting  these  equal  to  zero.   The  constraint 
equation  is  also  imposed. 

The  following  set  of  equations,  as  developed  in 
Appendix  B-4,  results.   The  equations  must  be  numbered  as 
shown  to  facilitate  enforcem.ent  of  boundary  conditions. 
This  is  described  later  in  the  development. 

:      h    .1^   ^IJ  ^j  -^oE^r    0  (1) 

4  D 

A.   E   k«   V.  -  P   f4  =  0  (2) 

1  j^]^   2j   J     o  EC 

4       _  i|     _ 

(llll)    j  =  i   3j   (j+i-4)    (i/2)  .^    Ij  (j+i-2) 

-  Q  ^^-^^  =  0  (i-1)   ' 

EC 

n  _  H       _ 

\k:l^    ,.fl  %•  ^(J  +  i-^)'  '(1/2)  jfi^^2J^(j.i-2) 

p  5. 

-  -f^  =0        i  =  4,6, ...2n  (1) 

A    Z   k^  V,.,,       -^liHill  =  0  (2n+l) 

n  .  -^   3j   (j+2n-2)       EC 
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n 

V 

E      A.    • 

-      °    =    0 

1=1      ^ 

A 

4    i|     _        _ 
Jl  .fl^lj^(j  +  2m-2)^(l+2r.-2)-  "EC  =  ^         (2n+2+m) 

m  =  1,2,.  . .  ,n 


(3n+3) 


where : 

-      th 

D.  =  i    component  of  the  global  consistent  load 

distribution  vector,  5D^„  ,_v 

'  lx(2n+2) 

2 
A   =  strength  constant  S  /6E  as  developed  in 

Section  II.         ^ 

V.  =  i    component  of  the  global  displacement 
vector. 

This  is  a  set  of  3n+3  nonlinear  algebraic  equations 

which,  when  solved  simultaneously,  results  in: 

i)  values  for  n+1  nodal  displacem.ents  and  n+1  nodal 

slopes  at  maxim.um  load  Intensity. 

ii)  the  n  optimum  element  areas 

iii)  the  maximum  load  intensity,  P 

No  boundary  conditions  are  included  in  the  preceding 

global  formulation.   Since  the  derivatives  may  only  be  taken 

with  respect  to  variables,  a  fixed  displacement  or  slope 

boundary  condition  requires  the  removal  of  the  equation 

associated  \-;ith  the  variation  of  T*  with  respect  to  that 

boundary  condition.   Its  value  must  be  substituted  in  the 

remaining  equations.   Boundary  conditions  must  be  consistent 

with  the  finite  elem.ent  model,  i.e.  they  m.ust  occur  at 

nodal  points . 
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For  example,  if  the  j    state  variable,  v.,  is  fixed, 
the  J  "  equation  is  eliminated  from,  the  set,  and  the  value 
for  V.  substituted  in  the  remaining  equations. 

In  this  way,  boundary  conditions  are  included  in  the 
formiulation  for  specific  problems. 

The  set  of  equations  with  boundary  conditions  considered 
is  the  finite  element  model  for  the  approximate  structural 
optimization  of  a  particular  rectangular  cross-section  beam 
of  uniform  height  and  fixed  volume. 
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IV.   APPLICATION  OF  THE  FINITE  ELEMENT  MODEL 

In  order  to  apply  the  model  developed  in  the  preceding 
section  we  must  find  some  method  for  solving;  the  system  of 
nonlinear  simultaneous  algebraic  equations.   The  size  of  the 
system  for  even  a  small  number  of  elements  requires  the  use 
of  a  computer  algorithm  for  solution.   The  only  such 
algorithm  available  to  the  author  was  a  Fortran  IV-G  sub- 
routine called  SUBROUTINE  NLNSYS .   This  subroutine  is 
included  in  the  computer  program  section  of  this  report. 

A.   SUBROUTINE  NLNSYS 

This  subroutine,  in  theory,  solves  simultaneously,  by 
an  iterative  scheme  developed  by  Brown  [8],  [9])    [10],  a 
system  of  nonlinear  algebraic  equations  of  any  order.   As 
it  is  dimensioned  it  is  restricted  to,  at  most,  100 
equations , 

NLNSYS  requires  as  input  a  vector  of  starting  values 
for  all  unknowns.   This  starting  guess  vector  must  be  in 
the  region  necessary  for  convergence  or  the  process  will 
diverge.   This  region  cannot  be  defined  and  different  start- 
ing guesses  may  be  necessary  to  obtain  convergence.   It 
should  be  noted  that  large  systems  of  equations  require 
large  starting  guess  vectors.   Since  the  probability  of 
a  component  being  outside  the  region  of  convergence  increases 
with  the  number  of  com.ponents,  the  probability  of  a  starting 
guess  vector  being  within  this  range  decreases  with 
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increasing:  system  size.   This  v;as  found  to  be  a  difficulty 
that  limited  the  size  of  the  example  problems  to  three 
elements . 

NLNSYS  also  requires  an  external  subroutine  that  speci- 
fied the  system  of  equations.   Because  the  system  of 
equations  developed  in  Section  III  changes  with  the  number 
of  elements  and  boundary  conditions  of  a  problem,  a  general 
subroutine  is  desired  that  will  generate  the  proper  set  of 
equations  for  a  given  problem.   SUBROUTINE  FCNLST ,  provided 
in  the  program  section,  accom.plishes  this  task. 

B.  SUBROUTINE  FCNLST 

SUBROUTINE  FCNLST  generate   a  set  of  equations  which  are 
evaluated  at  a  point  provided  by  NLNSYS.   NLNSYS  selects 
only  one  equation  at  a  time,  thus  the  set  of  equations  must 
be  ordered  in  some  manner. 

SUBROUTINE  FCNLST  sequences  the  equations  for  a  given 
number  of  elements  as  they  are  ordered  in  Section  III.   It 
evaluates  these  equations  at  the  point  specified  by  NLNSYS. 

The  equations  derived  from  parameters  that  are  fixed 
by  boundary  conditions  are  removed  from  the  general  set  and 
NLNSYS  operates  on  one  of  the  equations  in  this  reduced  set. 

Problem  parameters  are  read  into  the  program  in  SUBROUTINE 
FCNLST  the  first  time  that  it  is  called  by  NLNSYS. 

C.  SUBROUTINE  OPBEAM 

SUBROUTINE  OPBEAM,  as  provided  in  the  program  section, 
reads  the  initial  guess  vector  into  the  program  and  provides 
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the  output  from  the  program.   It  permits  an  efficient  means 
of  trying  multiple  starting  guess  vectors  for  a  given 
problem.   It  starts  the  solution  to  the  problem  by  calling 
NLNSYS . 

D.   INSTRUCTIONS  FOR  PROGRAM  USE 

It  is  first  necessary  to  define  the  problem  to  be 
solved.   We  specify: 

1)  the  length  of  the  beam  (in)  -  BLENG 

2)  the  volume  (in^)  -  VOL 

3)  the  height  (in)  -  HI 

^)   the  modulus  of  elasticity  (lb/in  )  -  E 
5)   the  yield  stress  (Ib/in^)  ~  SIGYP 
Selecting  the  number  of  elements,  NI ,  we  number  the  dis- 
placements and  slopes  at  element  boundaries.   This  must  be 
done  as  indicated  in  Figure  ^. 


„  =  slope  v^  =  slope   v^  =  sloDe   Vp=slope   '^'2'NI  +  2 


Figure  H.      The  Finite  Element  Numbering  Scheme 
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We  consider  the  boundary  conditions  of  the  problem  and 
form  two  vectors.   One,  called  IPC,  is  a  vector  of  subscripts 
of  the  components  of  the  displacement/slope  vector,  V,  which 
are  fixed.   These  are  ordered  in  increasing  value.   The 
second,  called  BC ,  is  a  vector  of  the  values  of  these  fixed 
components .   The  number  of  boundary  conditions  is  called 
NBC.   For  example,  if  the  displacem.ents  at  the  left  and 
right  ends  of  the  beam  are  fixed  at  zero  and  the  slopes  at 
these  points  are  variable  (simply  supported),  the  vectors 
are  : 

IBC(l)  =  1  BC(1)  =  0 

1BC(2)  =  (2  X  NI)  +  1      BC(2)  =  0 

The  total  number  of  parameters  is  equal  to  (3  x  NI )  +  3 
and  is  denoted  NT.   Subtracting  the  number  of  boundary  con- 
ditions, NBC,  v;e  obtain  the  total  number  of  unknowns, 
NU  =  NT  -  NBC.   This  is  also  the  size  of  the  reduced  system. 

From  the  given  loading  condition  on  the  beam  we  form  the 
global  consistent  load  distribution  vector  <D>  for  the  number 
of  elements  chosen.   This  is  denoted  CLOAD  in  the  program. 

The  initial  guess  vector,  X,  is  form.ed  by  guessing  the 
magnitude  of  the  unknown  displacements  and  slopes,  the 
areas  of  the  elements,  and  the  maximum  load  intensity. 
These  must  be  ordered  as  f  ollov/s  : 

1)  Displacements  and  slopes  as  ordered  in  the 
vector  {v}  v;ith  boundary  conditions  removed. 

2)  Element  areas  starting  at  the  origin. 

3)  Maximum  load  intensity. 
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V/e  select  values  for  MAXIT,  NUMSIG,  and  IPRIIIT  from 
their  descriptions  in  SUBROUTINE  NLNSYS. 

A  calling  program  is  v/ritten  v;hich  dim.ensions  X 
(the  starting  guess  vector).   This  dimension  is  NU .   The 
program  calls  SUBROUTINE  OPBEAM  with  numerical  values  for 
the  arguments  except  for  X.   A  sample  calling  program,  is: 
DIMENSION  X(NU) 

IL  U.  it  ii 

CALL  OPBEAM(NU,MAXIT, NUMSIG, IPRINT,X) 

STOP 

END 

* (numerical  values) 

The  user  must  provide  the  two  statements  for  SUBROUTINE 
FCNLST  as  described  in  the  program.   The  first  is  the 
dimension  statement: 

DIMENSION  IBC(NBC) ,BC(NBC) ,CLOAD( 2xNI+2 ) ,G(NT) , 
H(NT) ,V(2xNI+2) ,A(NI) 
The  dimensions  are  numerical  values. 
The  second  is  the  data  statement: 
DATA  NU/«/,NBC/«/,NI/«/,NT/«-/ 
where  *  indicates  the  numerical  values  for  the  parameters. 
Data  cards  are  prepared  and  inserted  in  the  following 
order  and  format : 

1)  The  guess  vector  X  -  6E12.5 

2)  The  physical  constants: 
BLENG,E,HI,VOL,SIGYP  -  5E12.5 

3)  The  boundary  conditions: 
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IBC(l)  ,  E'C(l)  -  I10,E12.5 
IBC(2) ,  BC(2)  -  I10,E12.5 

IBC(NBC) ,BC(NBC)  -  I10,E12.5 
^)   The  consistent  load  distribution  vector 
CLOAD  -  6E12.5 

The  final  output  may  be  one  of  the  following: 

1)  A  singularity  was  generated  at  the  output  vector 

2)  The  iteration  limit  was  reached. 

3)  A  solution  v/as  reached. 

In  case  1  a  new  starting  guess  should  be  miade  and  the 
problem  re-run.   In  case  2,  if  the  number  of  iterations  was 
small  (<30),  try  a  larger  number  of  iterations.   If  it  was 
not  small,  try  a  new  starting  guess.   Many  guesses  may  be 
required  before  a  solution  is  reached. 

An  example  problem  set-up  is  shown  in  Appendix  D  along 
with  the  results  obtained  for  different  problems. 
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V.   RESULTS,  CONCLUSIOTTS,  AND  RECOMMENDATIONS 

A.   RESULTS 

Figures  22,  26,  and  28  compare  the  three  element  finite 
element  optimum  shape  with  the  exact  optim.um  shape  for  a 
specific  beam  under  a  specific  loading  condition.   Table  I 
compares  the  maximum  load  intensities  for  the  finite  element 
designs  of  Figures  22-28  with  those  for  the  beam  of  equal 
volume,  length,  and  height  but  constant  cross-section. 

TABLE  I 
COMPARISON  OF  MAXIMUM  LOAD  RESULTS 


Loading 

A 

B 

C 

D 

■  ■ 
E 

Simply  Supported 
Uniform  Load 

510  lb/in 

764  lb /in 

721 

lb/in 

l4l 

5.6 

Concentrated 

Loads  at 

X  =  10  in 

100  KIP 

— 

194 

KIP 

194 

— 

X  =  20  in 

55  KIP 

— 

105 

KIP 

191 

- 

X  =  30  in 

40.7  KIP 

— 

76 

KIP 

187 

- 

X  =  40  in 

34.4  KIP 

— 

63 

KIP 

183 

- 

X  =  50  in 

31. 4  KIP 

— 

58 

KIP 

181 

- 

X  =  60  in 

30.6  KIP 

- 

57 

KIP 

186 

G.e 

Triangular  Load 

996  lb/in 

- 

1830 

KIP 

185 

- 

Cantilever  Beam 
Uniform  Load 

191  lb/in 

382  lb/in 

360 

lb/in 

188 

5.6 

A  =  Maximum  Load  for  Constant  Area  Beami 

B  =  Maximum  Load  for  Exact  Optimum  Beam 

C  =  Maximum  Load  for  Three  Element  Beam 

D  =  ^  Increase  of  Load  from  Case  A  to  Case  C 

E  =  ^  Deviation  of  Case  C  from  Case  B 
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It  is  clearly  evident  in  Figures  22,  26,  an:  1:  that  the 
three-elenent  design  approximates  the  exact  opti:-..-.  shape. 
The  best  criterion  for  the  accuracy  of  the  appr: "i~ation , 
however,  is  not  the  bean's  appearance.   Since  o-jt  ~oal  is  to 
maximize  the  load  intensity,  the  accuracy  of  th-  ipproxim.a- 
tion  is  best  measured  by  the  closeness  of  the  Icii  intensity 
allowed  by  the  finite  elem.ent  design  to  that  all:-:ed  by  the 
exact  optim^um  design. 

In  those  examples  where  exact  solutions  were  >nov;n,  even 
the  coarse  three-element  approximation  allows  a  ~axim.um  load 
intensity  that  is  close  to  that  of  the  exact  op*:i""um  shape. 
Table  I  shows  that  the  three-elem.ent  design  gives  maximum 
loads  within  5.6%  of  the  exact  design  maximum  f  c  r  the  sim.ply 
supported  and  cantilever  beams  under  uniform  lo?. :,  and  within 
6.6%  for  the  simply  supported  beam  under  a  c on cer.t rated 
load  at  its  center. 

Attempts  were  made  to  increase  the  number  of  elements 
in  the  simply  supported  beam  under  a  uniform  lc?,.r. .   This 
was  desired  in  order  to  show  that  the  shape  woul :;  converge 
to  the  exact  solution  with  an  increasing  number  of  elements. 
However,  difficulty  was  encountered  in  finding  an  initial 
guess  vector  within  the  range  of  convergence.   lue  to  the 
limited  tim.e  available  and  the  need  to  check  the  model's 
validity  for  other  loading  and  boundary  conditirr.s,  the 
three-element  m.odel  v;as  used  throughout. 
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B.   CONCLUSIONS 

The  above  results  show  that  the  finite  ele:---r.t  model 
developed  in  Section  III,  for  the  structural  07.  t  imization 
of  a  homogeneous  beam  of  uniform  fixed  height ,  =.:-.d  fixed 
volume,  is  valid.   Furtherm.ore ,  the  results  sh:--  that  a 
relatively  small  num.ber  of  elements  may  be  usei  to  design 
a  beam  that  will  closely  approximate  the  strer.eih  of  an 
exact  maximum  strength  design. 

As  a  practical  m.atter,  the  construction  of  a  beam,  with 
the  shape  of  the  finite  element  solution  is  miuz-h  easier 
than  constructing  the  exact  optimum  shape.   Ir.  the  exact 
optimumi  design,  points  v/here  the  cross-sectior.al  area  is 
zero  must  have  some  material  added  to  account  for  the  shear 
effects  that  were  neglected  in  the  formulation.   Much  less, 
if  any,  material  must  be  added  to  the  finite  element  design 
at  the  same  points.   The  finite  element  design  would  have 
stress  concentrations  at  the  elem^ent  boundaries  which  are 
not  present  in  the  exact  optimum  design. 

The  finite  element  model  provides  a  means  rf  approxim.at- 
Ing  solutions  to  optimization  problems  which  are  impossible 
to  solve  in  a  closed  form  by  classical  techniques  and 
extremely  difficult  to  solve  numerically.   This  is  due 
mainly  to  the  algebraic  nature  of  the  finite  element  method 
and  the  differential  equation  nature  of  the  classical 
methods . 

The  model  developed  was,  of  necessity,  fcr-  the  sim.plest 
type  of  beam  and  loading  conditions.   The  approach  used. 
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however,  should  be  valid,  v/ith  some  modifications,  for  the 
more  comiplex  cross-sections.   It  has  been  shown  that  the 
finite  elem.ent  techniques  can  be  successfully  applied  to 
the  problems  of  beam  optimization. 

C .   RECOMMENDATIONS 

Many  more  problems  for  v;hich  the  developed  model  applies 
should  be  solved  with  various  combinations  of  loading  and 
boundary  conditions  for  a  parametric  study  of  the  beam 
optimization  problem.   Problem.s  with  more  elements  should  be 
solved  to  show  convergence  to  known  exact  solutions. 

Although  SUBROUTINE  NLNSYS  provided  solutions,  the 
algorithm  is  very  sensitive  to  the  starting  guess  in  rela- 
tion to  the  region  of  convergence.   More  than  one  attempt  is 
usually  required  to  obtain  a  solution  to  a,  given  problem. 
An  algorithm  that  does  not  exhibit  this  sensitivity  is 
required  to  enable  m.ore  efficient  problem,  solving  using 
finite  element  techniques.   At  the  present  time  very  few 
algorithms  for  solving  systems  of  nonlinear  equations  exist. 

Another  technique  that  is  possible  would  be  to  use  some  - 
type  of  direct  search  or  steepest  descent  algorithm  to  find 
an  optimum  solution  to  the  finite  element  potential  energy 
functional  itself.   VJith  this  type  of  formulation  it  may 
be  possible  to  allow  variable  element  areas.   A  solution 
might  be  more  readily  attained  as  compared  with  a  solution 
to  the  system  of  nonlinear  equations.   Som.e  work  v;as  started 
in  this  direction  late  in  the  investigation.   The  initial 
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results  looked  promising'  but  tlr.e  restrictions  dictated 
that  this  approach  be  terminated  in  favor  of  com.pletinp 
the  presented  formulation. 

Models  for  other  cross  section  types  could  be  developed 
using  a  similar  approach  to  that  presented  in  Section  III. 
Also  models  xvith  other  than  equal  size  elem.ents  or  for 
beams  of  m.ore  than  one  material  would  be  of  interest. 

In  the  laboratory,  an  experim.ent  could  be  run  comparing 
the  behavior  of  an  exact  optimum  beam  with  that  of  its 
finite  element  approximations  to  see  if  the  results  compare 
with  those  predicted  by  the  models. 

The  entire  area  of  the  application  of  finite  element 
methods  to  optimization  has  hardly  been  touched.   The 
finite  element  m.ethod  is  a  powerful  tool  and  more  research 
should  be  done  in  its  application  in  this  field. 
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APPENDIX  A 

EXAMPLES  OF  BEAM  STRUCTURAL  OPTIMIZATION  PROBLEMS 
USING  THE  VARIATIONAL  FORMULATION 

1.   Consider  a  simply  supported  beam  of  rectangular  cross- 
section  with  constant  height,  h,  and  variable  width,  b(x). 
The  beam  is  under  a  uniform  load  distribution  and  has 
given  volume  and  length,  V   and  L  respectively.   Find  the 
optimum  shape,  b(x)  such  that  the  load  intensity  P   is  a 
maximum.   The  beam  has  a  modulus  of  elasticity  E  and  yieli 
strength  S  . 

y 

Solution : 
The  moment  of  inertia  of  the  cross-section  is: 

12  2 

T    bh     h   /Ku^    h   .  , 

^  =  12-  =  12  ^^^^  =  12  ^ 
Thus  C  =  h  /12  and  n  =  1.   Using  equation  (11-8)  with  n  =  1 


obtain 


V  (Lx-x^) 
A  =   ° 


L      2 
/  (Lx-x  )dx 


0 

Invoking  equation  (11-10)  to  obtain  (11-11)  the  value  of 
becomes : 

^  =  61 

Substituting  the  values  for  C,  X,  and  n  into  equation 
(11-9)  yields: 


kl 


hS  V 

P    =  y  o 


o     L 

3/   (Lx-x'')dx 
0 


The  integral  can  be  evaluated  in  closed  form 
/   (Lx-x'^)  dx  =  L^/6 


Thus  : 


0 


6V 

A(x)   ,  =  ~  (Lx-x^) 
opt   ^3 


v/hich    gives  : 


6V 


b(x)       ,     =    — ^    (Lx-x^) 
opt        ^^3 

The   beam   shape    is   graphically  -depicted   in  Figure    5.      The 
maximum   load   intensity    is: 

2hS   V 

P  =    ^ 

o  T  3 

max  L-^ 

This  shows  that,  for  the  given  problem,  maximum  load 
intensity  varies  linearly  with  the  height,  yield  stress, 
and  volume,  and  inversely  with  the  cube  of  length. 

2.   Consider  a  simply  supported  beam  of  rectangular  cross- 
section  with  a  constant  width,  b,  and  variable  height,  h(x) 

The  beam  is  under  a  uniform  load  intensity  P   and  has  a 

^   o 

given  volume,  V  ,  and  length,  L.  Find  the  optimum  shape, 
h(x)  such  that  the  load  intensity,  P  ,  is  a  maximum.  The 
beam  has  a  modulus  of  elasticity,  E,  and  yield  strength,  S 


h2 


Solution : 


The  mordent  of  inertia  of  the  cross-section  is 


I  = 


bh- 


12        2 


1   a3 


Thus  C  =  X-  and  n  =  3.   Using  equation  (II-8)  with  n 

12b 


=  3 


A  = 


o 

/  (Lx-x  )  dx 
0 


Equation  (11-13)  gives 


X  =  S  /2E 

y 


With  A,  C,  and  n  above,  equation  (II-9)  becomes 


?  V 

P     =    y   o  . 

o      3b 


1 


r  L 


2  '-2 
/  (Lx-x  )  dx 


L  0 


The  integral  may  be  evaluated  and  is 


/   (Lx-x  )  dx  =  — Q— 


0 


Thus 


A 


8V  (Lx-x^)      2.55  V  (Lx-x^) 
o  o 


opt 


ttL' 


2 


Solving  for  h(x)  : 


h(x)  0. 


2  '^ 
2.55  V^(Lx-x^) 


bL 


2 
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h(x)  Is  gra|±ilcally  depicted  in  Figure  6.   The  maximum  per- 
missible load  intensity  is: 

64s  V^    2. ITS  V^ 

Po    =  ^L^   =.  1-^ 

max     2,  T  l\  ,  -.  4 

s   bL      bL 

Here  Po    varies  linearly  with  yield  stress,  para- 
bolically  with  volume  and  inversely  with  width  and  the 
fourth  pov/er  of  length. 

3.   Consider  a  simply  supported  beam  of  circular  cross- 
section  under  a  uniform  load  distribution.   It  has  a  given 
volume  and  length,  V   and  L  respectively.   Find  the  optimum 
radius,  r(x)  such  that  the  load  intensity  is  a  maxim.um. 
The  beam  has  a  modulus  of  elasticity,  E,  and  yield  strength. 


S  . 

y 


Solution 


The  moment  of  inertia  of  a  circular  cross-section  is: 
4  2 

Thus  C  =  T—  and  n  =  2.   V/ith  n  =  2,  equation  (11-8)  becomes 

P  2/3 

V  (Lx-x  ) 
A  =    ° 


L      2  2/3 
/  (Lx-x  )    dx 


•0 

Equation  (11-12)  gives  the  value  of  A  for  a  circular  cross- 
section  as  : 

S  2 
A  =   ^ 


i|E 
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With  A,  C,  and  n  as  above,  equation  (3-9)  becor.es 


P   = 
o 


S  V 

y  o 


2/3 


1 


2  A 


r  L    2  2/3- 

/  (Lx-x  )    dx 

Lo 


3/2 


The  integral  may  be  evaluated  in  closed  form.   However  it 
involves  the  Gamma  Function  which  must  be  approximated. 
The  integral  is,  then: 

L     p  2/3  - 

/  (Lx-x"")    dx  ^  .0979  *  (L)  '^^ 
0 


Thus 


A 


V  (Lx-x^) 
o       ' 


2/3 


opt  -   .0979l'^/3 


2/3 


V^(Lx-x^) 


308  L 


7/3 


The  beam  shape  is  depicted  in  Figure  7-   The  maxim.um 
load  intensity  is  given  by: 


Po 


S  V 

y  o 


3/2 


max 


.1067  L 


7/2 


Po    is  linear  in  yield  strength,  proportional  to  the 
max  "^  to   >  i-   t- 

7/2  power  of  volume  and  inversely  proportional  to  the  7/2 
power  of  length. 

Reference  [11]  gives  optimum  beam,  shapes  for  example 
problems  1  and  2  derived  in  a  different  manner.   The  shapes 


^5 


shown  in  Figures  5  and  6  are  identical  to  those  shapes 
The  maximum  load  intensities  also  correspond. 


H 


y 


i i. 

h 


± 


L 


♦-  X 


rfyrr 


-^ 


?v^ 


i.no 


p(x)  =  + 


(Lx-x  ) 


«•  •nr^^b 


ODtlrur  Pear  Sha.De  for  a  Pectang-ular  Cross- 
section  Peam  of  Uniforr.  Height  unrier  a 
Uniform  Load. 


n 


uniform  load  distribution 


intensity  P 


y 

Figure  6 . 


Optirrum  Peam  Shape  for  a  Pectanp-ular  Pross- 
section  Pear  of  Uniforr  V'idth  under 
Uniform  Load. 
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1.50 


-.75 


-1.00  — 


-1.50 


uniform  load  distribution 


intensity  ^ 


,:       t»  X 


o 


y 


Figure  7.   Optirur  ^esr.  Sbane  for  a  Circular  ProsS' 
section  Pear  under  a  Uniform  Load. 
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APPENDIX  P 


CALCULATIONS  FOR  THE  FINITE  ELET^ENT  I^^THOD 
FOR  PEAF  STRUCTURAL  OPTIMIZATION 


1 .   The  Cubic  Displacement  Function 

VJe  wish  to  find  a  cubic  displacement  function  vector 


<N(x.  )>  such  that 
1 


v(x^)  =  <N(x^)>  J  .^ 


where  v(x,)  Is  a  cubic  polynomial. 


2^1 


Figure  8.   The  Deflected  Element. 


Q      2 
We  let  v(x.)  =  ax.  +  6x,  +  yx..  +  6  over  the  Interval 

0  <  X.  £  A.   Boundary  conditions  are: 
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(1)  v(0)  =  v^ 

(2)  v'(o)  =  e^ 

(3)  v(A)  =  v^ 
(^)  V' (A)  =  02 

Substituting  boundary  condition  (1)  in  v(x.)  gives  6  =  v. 
Taking  the  first  derivative  of  v(x.)  with  respect  to 


X.  yields : 


v'  (x. )  =  3ax?  +  26x.  +  y 
1        1       1    ' 

Boundary  condition  (2)  in  v'(x.)  yields: 

Substituting  the  values  for  y  and  5  into  v,  and  v'  and 
invoking  boundary  conditions  (3)  and  (4)  yields: 

2v^   01   2V2   e^ 

«  =  —-  +  —  --—+  — 

A^     A     A^     A 

3v    20     3v    e 
^  "  "  a2   -   A   ^  .2   -  A 


A' 
In  vector  notation: 


/  2   1    2   1  X  /  \ 


■a3  a2    a3  a2 


<- i  -  i  j^  -  iX'h 


y    =    <0    1    0    0>    {v}. 
6=<1000>{v}. 


Thus  in  vector  notation 
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v(x.) 


2x; 


3x: 

2 


2x' 


+  1 


,  2 


+  X 


2x: 


^4\ 


i  2\ 

X.    X.  K 


Therefore  the  shape  function  vector  <N(x.)>  is: 


2x 


<N(x. )> 
1 

Ix^ 


1 


a3 

3 
i_ 
2 


3x?    \/x?    2x? 
2 


2x 


3x: 


The  second  derivative  of  this  vector  v/ith  respect  to  x.  is 

'^         1 


12x 


<N(x. )  > 
1 


r     \  ; 6x. 
6  W   1 


a2   a2 


12x 


i+  f 


6x .    ^ 
1    2 


a2  U2 


2.   The  Modified  Element  Stiffness  Matrix 

The  modified  element  stiffness  matrix  is  defined  as 


A 


T 


[k«]  =  /   <N"(x. )>   <N"(x. )>  dx. 

1  1      1 

4x4     °    4x1       1x4 


Taking  the  indicated  vector  product  and  integrating  yields 


M  = 


12/A 
6/A' 
-12/A 
6/A' 


3 


6/A 

4/A 

•6/A 

2/A 


2 


-12/A 

-  6/A' 
12/A 

-  6/A 


3 


3 


2 


6/A 

2/A 

■6/A 

4/A 


2 


3.   The  Consistent  Load  Distribution  Vector 


th 


Assume  that  the  load  per  unit  length  over  the  i    element 
may  be  written  as  an  r   order  polynomial 
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_    -      _  p  -  r 

P(x.)  =  a  +  3x.  +  yx.  +  ...  +  nx 

where  a,  6,  Y,...,ri  are  constants  v;hich  are  deternined  from 

known  values  of  load  intensity  at  r+1  loactions  in 

0  <_  X.  £  A.   These  known  values  are  represented  as  fractions 

of  the  maximum  load  intensity,  P  ,  over  the  beam  (n<x<L) 

P(Xi) 

such    as    a   =    a    •    — ^=^ where    0    <    x.    <    A.      Thus: 

Po  -     1    - 

P(x.)    =    P    [a   +    Bx.    +    yx?    +    ...    +    r]x^] 
1  o  1  '     1 

The    load    on    an    increment    dx.    of   the    i         elem.ent    is: 

1 

dP(x.)    =    P    [a   +    Bx.    +   yxT    +    ...    +    nx^]    dx. 
1  o"-  1  'l  1  1 

The   work    done    on   the    i         element    by    this    increm.ental    force 
is: 

dw.    =    dP(x. )    v(x. ) 
1  1  1 

but : 

v(x. )    =    <N(x. )>    {v} . 
1  11 

Ix^  4x1 

Thus  : 

?  r 

dw.    =    P    [a+Bx.+yx.  +  .  .  .+T1X.  ]    <N(x.)>    {v}^    dx^ 

th 
The   work   done    on    the    i         element    is: 

A 

w.    =    /      dw. 
^         0 


w.    =    /      P    [a+Bx.+yx.+. . .+nx. ]    <N(x.)>    (v) .    dx. 
1  „  oil  1  1  1  1 


53  • 


Letting  [a+Bx^  +  Yx^.  +  .  .  .  +  nxr]  be  denoted  by  0,  ;;,;  and 
noting  that  P^  and  {v}^  are  independent  of  x.  ,  t'r  ■=:  above 
expression  becomes: 

A 

w.  =  P   /   0  (x. )  <N(x. )>  dx.  {v}. 
1     o      ri       i'     11 

The  integral  is  defined  as  the  element  consi^Tent  load 
distribution  vector  <D>. 


<D>.  =  /   0  (x. )  <N(x. )>  dx. 
1    rs        r'l'   ^,1      1 

1x4     0  1^^' 

The  work  done  on  the  i    element  may  then  be  wriTTen 

w.  =  P   <D>.  {v}. 

1      O      11 

The  work  done  on  the  beam  is 


n  n 

W  =  1      w.  =  P    E   <D>.  {v} . 

.   T    1       o   .   T       11 

1=1         1=1 

It  may  be  noted  there  is  no  requirement  for  the  -lement 
load  distribution  polynomial,  P(x.)j  to  be  the  s  ?.~"e  for  each 
elem.ent.   On  the  contrary,  each  element  may  have  a  different 
load  distribution  polynomial,  P(x.),  and  the  disTribution 
may  even  be  discontinuous  at  element  boundaries.   The  only 
requirement  is  that  the  ratio  of  P(local)/P   be  '.-".r.own  at 
r+1  points  over  the  element  in  order  to  define  tr.f  r    order 
distribution  polynomial. 

We  take  the  follov/ing  case  as  an  example: 

Consider  the  first  three  elem.ents  of  an  r.  element 
beam.   Given  that  the  distributions  of  load  on  elements  1, 
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2,  and  3  are  cubic,  linear  and  quadratic  respectively,  and 
also  given  the  values  of  P  /P   where  z  =  A,  B,  C H, 

Z      O  3  i  i  3    3 

we  seek  an  expression  for  the  work  done  on  the  first  three 
elements  of  the  beam.   Figure  9  shows  a  graphical  represen- 
tation of  the  problem. 


4 


y  3V 


Figure  9-   Multiple  Loading  Types  on  a  Peam 


The  values  P„/P   through  Py^/P   define  the  cubic  distri- 
A   o  Do 

bution  function  O^Cx,  )  over  the  first  elem^ent.   The  values 

P,^/P   and  P-r-i/P   define  the  linear  distribution  function 
E   o       F   o 

Q-,(x„)  over  the  second  element.   The  values  of  Pt-,/P 
Id  F   o 

through  P  /P   define  the  quadratic  distribution  function 
HO  ^ 

Op(x^)  over  the  third  elem.ent.   These  functions  in  turn  de- 


fine consistent  load  distribution  vectors  <D-,  >  3  <E)p> , 


and 


<D^>  for  the  elem.ents 


<D  >  =  /  0  (x. )  <N(Xn )>  dXn 

J-     Q  ^   ± 

A 

<T)^>    =  /  Q^(x2)  <N(Xp)>  dx^ 


<D  >  =  /   O^Cx^)  <N(x^)>  dx^ 
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The  work  done  on  the  first  three  elements  is: 

W(l-3)  =  Po[<D>i  ^v^}  +  <D>2  {v^}  +  <D>^  {v^}] 

M .   The  Concentrated  Load 

If  the  load  is  concentrated  vice  distributed,  the 
consistent  load  distribution  vector  is  formed  in  a  manner 
parallel  to  that  of  the  distributed  load. 

We  apply  a  concentrated  load  of  intensity  gP  C0<  6<_1) 
at  a  distance  aA(0  ±  ol   <_  1)    from  the  origin  of  the  i"*^^ 
element . 


>•  X 


Figure  10.   A  General  Element  with  Concentrated  Load 


The  work  done  on  this  element  is 


w.  =  6P  v(aA) 
1      o 


But 


v(aA)  =  <N(aA)>  {v} 
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Therefore : 

w.  =  6P   <N(aA)>  {v}. 
i      o  1 

If  we  define  the  eler^.ent  consistent  load  distribution 
vector  as : 

<D>.  =  3  <NCaA)> 
1 

th 
we  obtain  for  the  work  done  on  the  i    element  by  the  load 

P  : 
o 

w.  =  P   <D>.  {v}. 

10      11 

Note  that  the  expression  is  the  same  as  that  for  the  dis- 
tributed load.   However,  the  units  for  P   and  <D> .  are 

'  o         1 

different  -  P   is  in  units  of  force  vice  force  per  unit 
o  ^ 

length. 

If  multiple  concentrated  loads  are  applied  to  one 
element,  the  elem.ent  consistent  load  vector  is  the  sum  of 
those  obtained  by  considering  each  load  separately.   In 
symbolic  form: 

The  individual  load  vectors  are: 

<D^>^  =  B-^  <N(a^A)> 

<D.>.  =  B.  <N(a.A)> 
J  1     J      J 

The  element  consistent  load  vector  is: 

<D>.  =  E  e.  <N(a.A)> 
1    j   J      J 
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5 .   The  Finite  Element  Variational  Equations. 

V/e  have  as  the  global  augr.ented  potential  energy 
functional: 

T*  =  ^   E   A.  [v.}'^    [k*]  {v.}  -  P^   E   <D>,  {v.} 
S    ^  i=l   1    1  1      o  ^^-^     1    1 

n 
-  A   E   AA. 
1  =  1    ^ 

where  n  is  the  number  of  elements. 

Define  the  following  as  the  global  consistent  load 

distribution  vector   <D> 

lx(2n+2) 

12^1^2       H  2^^ 

<D>  =  <d;,  d,,  d,^  +  d:;-,  d,  +  d^,...,d  ,  +  d  ,  d- ,  d  > 

1'   1'   1     2'   1     2'    '  n-1    n'   n'   n 

1'  2'       3       2n+l'   2n+2 

<D>  is  formed  by  assembling  the  element  distribution  vectors 
as  indicated.   Thus  the  global  work  term  becomes: 

n        _ 

P    E   <D.>  (v. }  =  P   <D>         (v) 
o  .  _-,     1     1       o 

^  ^  lx(2n+2)  (2n+2)xl 

The  augmented  global  potential  energy  functional  that  will 
be  operated  on  is : 

T^  =  ^   E   A.  {v.}'^  [k«]  {v.}  -  P   <D>  {v}  -  XA   E   A, 

Taking  the  partial  derivative  of  T*  with  respect  to  each 
component  v.  in  {v}  and  setting  result  equal  to  zero  gives: 

J 

TT-^  =  0:   EGA,   E   k^.  v.  -  P  D,  =  0 
3v^  1  j=i   IJ   J     0  1 
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9T«  k  _  _ 

7~  =    0:   EGA,   Z   k«   V.  -  P  D„  =  0 
Sv^  1  j^^   2j   J     o  2 


8T* 


=  0:   EGA  .  ^    Z   k*   V, . 


+  EGA,./^^   Z   k*   v^.,.  TN  -  P  D,.  .N  =  0 
(i/2)  .  ,   Ij   (j+i-2)     o  (i-l) 


1  =  4,6,8,. . .2n 


8T«  4 

r-;^  =  0:   EGA  .  „    Z   k|t  .  V.  ...  ,,  ^ 


+  EGA, .  /^^   Z   k* .  V, .^.  „v  -  P  D. 
(i/2)       2j   (j+i-2)     o  1 


1  =  4,6,8,. . .2n 


8T*  4 

^ —  =  0:   EGA    Z   kn.  V, 


(2n+l)  j=l   -'^    '-  o   (2n+l) 

^:; ^- —  =  0:   EGA    Z   kf- .  v,.,„   _v   ^   -         „ 

This  gives  2n+2  equations. 

We  next  take  the  partial  derivatives  of  T*  with  respect 

to  the  A.'s  and  set  them  equal  to  zero: 
1  ^ 

-^  =  0:    Z    Z   k«-.  V         --  "^ 


V 


9A^    "•   j^^  ^:^  "ij  Mk+21-2)  ^(j+2i-2)   EG 

i  =  1,2,...  ,n 

This  gives  n  equations. 

At  this  stage  we  have  2n+2  equations  in  3n+3  unknov/ns 
hence  an  additional  equation  is  necessary.  This  equation 
is  the  volume  constraint : 
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n 
A  E   A.  -  V   =  0 

1  =  1   ^ 


Therefore  we  have  a  total  of  3n+3  equations  in  the  3n+3 
unknowns : 

2n+2  components  of  the  vector  {v} 

n  areas  A. 
1 

1  maximum  load  intensity  P 
V/e  therefore  have  a  mathematically  tractable  problem. 
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APPENDIX  C 
EXAMPLES  OF  CONSISTENT  LOAD  DISTRIBUTION  VECTORS 


1.   As  our  first  example  v;e  consider  a  three  element  problem 
under  a  uniform  load  distribution. 

P 


u4 


'                    1 

1 

p                   1 

'                v°                ' 

'               < 

1 

2 

? 

-  A H-» 


-•-  X 


^ 


Figure  11,   A  Three  Element  Beam,  Uniform  Load. 

Observing  the  general  elem.ent  of  Figure  12,  the  work  done 
on  this  element  is: 

A 
w.  =  /   P   <N(x. )>  dx.  {v} . 

1     ^  O       1       1      1 


y .  .  V 


Figure  12.   A  General  Element,  Uniform  Load 
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which  may  be  written: 

A 
w.  =  P   /   <N(x. )>  dx.  {v}. 

1      Op         1       1      1 

From  our  previous  definition  in  Appendix  B,  the  elem.ent 
consistent  load  distribution  vector  is: 

A 
<D> .  =  /   <N(x. )>  dx. 
^    0 

Evaluating  the  integral  v/here  <N(x.)>  is  as  defined  in 
Appendix  B: 


/  <M(x^)>  dx^  =  <|,  fi.  |,  -f!> 


Thus 


<D>   =  <^  ^  ^  _  ^  >  (C-1) 

Applying  this  to  the  three  elem.ent  problem  and  noting 
that  the  load  intensity  is  the  same  for  each  element,  we 
obtain: 

<D>^  =  <D>2  =  <D>3  =<|,  1^,  |,  -  I2  > 

We  assemble  these  element  vectors  to  obtain  the  global  con- 
sistent load  distribution  vector  <D> 

2  2 

<D>  =<|,  I2'  ^>    0'  A,  0,  |,  -  ^  >  (C-2) 

The  work  done  on  the  beam  is : 
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o   2 ' 


A   A 


A     A' 


V, 


j^.  A,  0,  A,  0,  2>  -  Y2  \  V 


V 


V 


6 


_7 
/8, 


2.   As  a  second  example  let  us  take  the  following  problem 


P  /2 
o 


>                                          ' 

P 

'  o 

" 

f                                           1 

] 

9 

Q 

P  /2 
o 


-t"Y*— — — '  A 


«- 


>»> 


y  ,v 


Figure  13-   Three  Element  Beam,  Multiple  Uniform  Loads. 

The  element  consistent  load  distribution  vector  for  a  uni- 
form load  from  equation  (C-1)  is: 

<D>   =  <  ^   ^   ^   -  ^  > 
^  i      2'  12'  2'    12 

In  this  problem,  hov/ever,  the  load  intensity  is  not  the 
same  on  each  element.   Therefore  the  element  consistent  load 
distribution  vectors  are: 
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Or: 


<D>^  =  <I^>,  =  2  <D>i 


<D>^  =  <D>. 
2       1 


<D>-,  =  <D> 


A   A^   A 


2 


3      ¥'  24'  l\'         2l\ 


<D>^  = 


<  ^  ^  _ 
2'  12'  2' 


A     A' 


12 


Assembling  these  vectors  we  obtain  for  the  global  consistent 
load  distribution  vector: 

-      A  A^   3^  A^  3A  :::A^  A  -a^ 
IT'  2?'  T~'  12'  T~'  12'  ^'  2T" 

3.   Let  us  now  look  at  the  following  problem  with  a  concen- 
trated load: 


Ub  «  ^^ 


o 


1 

p 

Q 

■♦«- 


-pisJ- 


y  ,v 


-^ 


Figure  l4.   Three  Element  Beam,  One  Concentrated  Load 
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Taking  a  general  element  with  a  concentrated  load 


Figure  15.   General  Element,  Concentrated  Load 


The  consistent  load  distribution  vector  for  this  element  is 


But 


Thus 


<D>.  =  <N(aA)> 


<N(aA)>  =  <(2a^-3a^+l) ,  (a^A-2a^a+aA) ,  (-2a^+3a^), 
(a^A-a^A)  > 


<D>,  =  <(2a3-3a^+l),  [ (aA ) (a^-2a+l) ] ,  [ (a^) (-2a+3) ] , 


[(a^A)(a-l)]  > 


(C-M 


In  the  specified  problem  there  are  no  loads  on  elements  2 
or  3  thus : 

<D>2  =  <D^>    =  <0,  0,  0,  0> 

The  elem.ent  consistent  load  vector  for  element  1  is  there- 
fore : 


<n>   -  <nr3A^>   -  ^  5    3A   27   9A  . 


(C-Ha) 
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Assembling  '^D>-,,  <D>2',  and  <D^>  we  obtain  the  global  con^ 
sistent  load  distribution  vector  <D>. 


<D>  =  <  ^2'  riT'  §2'  -  ^'  °'  0,  0,  0  > 


(C-5) 


h.      Let  us  take  example  3  and  add  an  additional  load  on 


element  1  of  intensity  P   a  distance  jr   from  the  left  e 


nd 


A 
4, 

P 

0 

P 
0 

r 

1 

2 

■3 

}*» A ^ 


-H^ 


Figure  16. 


Three  Element  Beam,  Tv;o  Concentrated 
Loads  on  Element  One. 


The  element  consistent  load  distribution  vector  for  a 

concentrated  load  P   a  distance  A/4  from  the  origin  is: 

o 


<D(^)>. 


=   <    il     9A   5_  _  3A  . 
32'  ^'  32'    6T 


(C-6) 


In  this  example  <D>p  and  <D>^  are  the  same  as  in  example  3 
That  is  <D>p  =  <D>^  =  <0,0,0,0>.  The  consistent  load  dis- 
tribution vector  for  element  1  is  the  sum: 


<D>^  =  <D(^)>  +  <D(^)> 


IT^ 


■IT' 
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Thus: 


Assembling  the  element  vectors  v;e  obtain  for  the  global 
vector : 


<D>  =  <1,  ||,  1,  -  ^,  0,0,0,0> 


(C-7) 


5.   For  this  example  let  us  take  the  following  beam  with 
concentrated  loads  as  shown: 


y  ,v 


Figure  17. 


Three  Element  Beam,  Concentrated 
Loads  on  Elements  one  and  Three. 


From  equation  (C-6)  we  have  for  <D> 


1 


<n>   -  <  27   aA   5_    3A  ^ 
^1  "    32'  6V  32'  "  FiT 


Element  2  has  no  loads  thus 


<D>2  =  <0,  0,  0,  0> 


Equation  (C-4a)  gives  the  element  consistent  load  vector 
for  a  loading  such  as  is  on  elem.ent  3  but  of  Intensity  P 
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The  element  consistent  load  distribution  vector  for  element 
3  is  therefore  one-half  of  that  value  in  (C-4a) .   Thus 


<D>^  =  < 


5    3A   27 


QA 


3     6T'  12H'  F^'  "  12^ 

Assembling  these  to  obtain  the  global  consistent  load 
distribution  vector: 


<5>  =  <  27   9A   5_    3A   5    3A    27     9A 
32'  PT'  32'  U\'  F^>  T2E>  W^'  ~  T2E 


(C-8) 


6.   As  a  final  example  let  us  take  a  three  element  beam 
under  a  triangular  load  as  shown: 


■S»x 


I 


A  -*«f^ A 


A  — — -*4 


y,v 


Figure  l8.   Three  Element  Beam,  Triangular  Load. 

We  first  examine  a  general  element  with  a  basic 
triangular  load  of  maximum  intensity  P  as  shown: 
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-fe-  X 


Figure  19.   General  Element,  Triangular  Load 


This  distribution  may  be  written 

Px. 
1 


P(x. )  =  —- 

1      A 


The  element  consistent  load  distribution  vector  by 
definition  is : 


A  X. 

<D>.  =  /   -v^  <M(x.  )>  dx. 
1    Q   A      1      1 


Evaluating  the  integral  we  obtain  for  the  basic  element 
consistent  load  distribution  for  a  triangular  load. 


P(x. )  =  Px. 
1      1 


<T)>   =  <  lA   A^   21A     Af_ 
i      20'  30'  FO"  '  ~  20 


(C-9) 


If  we  depict  the  exam.ple  problem  as  follows 


Figure  20. _   Three  Element  Beam,  Superposition  of 
Uniform  and  Triangular  Loads  .. 


y  jv 
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We  see  that  element  1  has  a  basic  triangular  load  distri- 
bution of  intensity  P  /3.   Element  2  has  a  uniform  load 
distribution  of  intensity  P  /3  with  the  triangular  distri- 
bution of  element  1  superposed.   Elem.ent  3  similarly  has 
this  superposition  but  the  uniform  load  intensity  is 

2Py3. 

From.  (C-9)  and  the  relation  P  =  P  /3  we  obtain: 

<D>   =  I  <D>^'^^^"S  =  <  A_  Af_   7A   _  A^  ^ 
^1    3    i  20'  90'  60'    20 

Prom,  the  result  for  <D>-,  above,  and  from  (C-1)  with 
intensity  P  /3,  the  uniform  portion  of  <D>p  is: 

2        2 

^uniform  _  1_  .  A   A_  A   -A 

u  2  3    2'  12'  2'  12 


Thus 


2        2 

^^^uniform  _  ^  A   A    A   -A   ^ 

^^"^2  =  ^  ^'  3^'  ^'  12~  "^ 


And 


<D>2  =  <D>^^i^"S  +  ^P^uniform 

_    13A   7A^   17A     2a2 
^^^2  "  ^  M~'  IHO'  ^0~'  -  115" 

The  third  element  consistent  load  distribution  vector  is 
developed  in  a  similar  m.anner: 

2         2 

^uniform  _2_A   A_A     A 

u  3  3    2"'  12'  2'  ~  12 

2         2 

_    A   A    A     A 

3'  1^'  3'    iH 


And: 


<-p,>   _  <.-p,v, uniform    .p.  triangular 
3      3  3 
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Or 


3     FcT"'  15 '  20'  ~  Wn 


Assembling  <D>-.  ,  <D>p,  and  <D>o  we  obtain 


<D>  = 


A_ 
20' 


a! 

QO' 


A   A 
3'  '^'■ 


2 

2A_   A 

3  '  ^ 


QA 
20' 


13A' 


>    (C-10) 
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APPENDIX  D 
FINITE  ELEMENT  OPTIFIZATIOM  EXAMPLES 

All  of  the  follov/ing"  examoles  consider  the  approximate 
structural  optiniization  of  a  rectangular  cross-section  beam 
of  uniform  height  and  fixed  volum^e.   The  physical  parameters 
are  : 

a)  beam  length  =  120  in 

b)  beam  height  =  6  in 

o 

c)  beam  volume  =  2200  in 

d)  modulus  of  elasticity  =  30x10   lb/in 

o       2 

e)  yield  strength  =  50x10   lb/in" 

A  three  element  approximation  is  used  throughout. 

1 .   Representative  problem  set-up 

As  a  representative  problem  consider  a  simply  sup- 
ported beam  under  a  uniform  load.  The  beam  is  depicted  in 
Figure  21.   NI  =  3. 


6  in 


"  '^  ''      • 

i         a         5 


/y  Ir/ 


v; 


1 


^0   in- 


<- 


-l\[)    in — -j^"? 


3 


•^0  in 


7 


2  'i4  ^^  '8 

Figure  21.   Three  Element  Problem,  Uniform  Load 
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Since  the  beam  is  sirrply  supported  the  boundary  condi- 
tions are  : 

V.  =  0 

V   =  0 

Therefore : 

IBC(l)  =  1      PC(1)  =  0.0 
IBC(2)  =  7      PC(2)  =  0.0 
NBC  =  2 

The  total  number  of  parameters  are: 

NT  =3x3+3=  12 

The  total  number  of  unknowns  are: 

'      NU  =  12  -  2  =  10 

From  Appendix  C,  the  consistent  load  distribution  vector  for 
a  uniform  load  on  three  elem.ents  is: 

2  2 

CLOAD  =  <  |,  1^, A, 0,^,0,    |,  -  ^  > 

Since  A  =  L/3  =  ^0  in,  CLOAD  becom.es: 

<  20,  133.333,  40,  0,  4o,  0,  20,  -133-33  > 

The  com.ponents  of  the  x  vector  are: 

X  =   <  v^,  v^,  v^,  v^,  Vg,  Vg,  A-^,  A^,  A^,  P^  > 


We  guess  values  for  the  x  vector 


^2 
^3 

= 

.03 
.9    in 

^2 

=    15 
=    28 

.    2 

m 

.     2 

m 

2 

^1) 

= 

.01 

A3 

=    15 

in 

^5 

= 

.9    in 

^6 

~ 

-..01 

^8 

—    —  • 

03 
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P   =  760  lb/in 
o 


VJe  select  values  for  f-^AXIT,  NUNSIG,  and  IPPINT.   Let  us 
require  50  iterations  maxinurn. 

MAXIT  =  50 
V/e  desire  four  significant  digits. 

NUMSIG  =  k 
We  desire  the  x  value  printed  at  each  iteration. 

IPRINT  =  0 

We  now  can  utilize  the  program.   The  calling  program  is: 

DIMENSION  X(10) 

CALL  OPBEA]VI(10,50,4,0,X) 

STOP 

END 

We  insert  the  following  cards  in  their  proper  locations  in 
FCNLST. 

DIMENSION  IBC(2) ,PC(2) ,CL0AD(8) ,G ( 12 ) ,H( 12 )  , 

V(8),A(3) 
DATA  NU/10/,NPC/2/,NI/3/, NT/12/ 

The  first  data  card  is  the  first  six  com.ponents  of  the 
X  vector  in  E12.5  format.   The  second  is  the  last  four 
components  of  the  x  vector  in  E12.5  f  orm.at .   The  third  card 
is  the  list  of  physical  parameters  BLENG ,E ,HI ,VOL ,SIGYP  in 
E12.5  format.   The  fourth  card  is  the  first  component  of 
IBC  and  the  corresponding  component  of  BC  (I10.,E12 . 5 )  • 
The  fifth  card  is  the  same  as  the  fourth  but  for  the  second 
components  of  the  vectors. 


7^ 


The  last  two  cards  contain  the  eig^ht  components  of  the 
consistent  load  distribution  vector  in  6E12.5  fornat.   The 
data  deck  is  thus: 

O.lE-1  0.9E0  -O.lE-1  -0.3E-1 

1.5E1  7.6E2 

6.0E0  2.2E3  5.0E4 


0.3E-1 

0.9E0 

1.5E1 

2.8E1 

1.2E1 

3.0E7 

1 

O.OEO 

2 

O.OEO 

2.0E1 

1.33333E2 

2.0E1 

-1.33333E2 

4.0E1  O.OEO  4.0E1  O.OEO 


This  program  was  run  and  the  following  solution  was 
returned  after  four  iterations. 

v„  =  .0310         A,  =  15.01  in^      P   =  720.6  lb/in 
2  1  P       o 

v^  =  .9566  in     Ap  =  2^.98  in 

v^  =  .0111        A^  =  15.01  In^ 

Vf-  =  .9566  in 
5 

Vg  =  -.0111 
Vg  =  -.0310 

We  are  interested  in  m.ainly  the  areas  and  maximum  load 
results.   If  the  areas  are  divided  by  the  height  (6  in)  we 
obtain  the  width.   Since  the  beam  is  symmetric  about  its 
centerline  we  may  plot  one  half  of  the  beam  to  depict  its 
shape.   Figure  22  is  a  plot  of  the  three  finite  element 
approximation  with  the  exact  solution  superim.posed . 

2  .   Other  exam.ple  problems 

The  problems  of  a  sim.ply  supported  beam  under  con- 
centrated loads  at  various  locations  were  also  solved  v;ith 
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the  computer  program.   The  results  for  these  problem.s  are 
shovm  in  Figures  23,  24,  25,  and  26.   Corresponding  loads 
to  the  right  of  the  beam,  m.idpoint  were  checked  and  the  beam, 
shapes  obtained  were  the  inverse  of  those  shown,  i.e.,  the 
half-v7idth  values  for  the  end  elements  were  interchanged. 

The  solution  of  a  simply  supported  beam  under  a  tri- 
angular load  distribution  is  shown  in  Figure  27. 

Finally  the  problem  of  a  cantilever  beam  under  a  uniform 
load  was  solved.   This  is  shown  in  Figure  28. 

Optimum  beam  shapes  for  the  simply  supported  beams  under 
1)  a  uniform  load  and  2)  a  concentrated  load  at  the  center 
and  also  for  the  cantilever  beam  under  a  uniform  load  are 
shown  in  [11]  .   These  exact  optim.um  shapes  are  superim.posed 
for  comparison  in  Figures  22,  -26,  and  28. 
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MAXIMUM  LOAD,  EXACT  SOLUTION 76^  lb/in 
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Figure  22.   ""hree  Elerent  ODtirum  Desip-n  for  a 
Simply  supported  Pearr  under 
Uniform  Load. 
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Loads  at  points  1  through  f   apply  in  Figures  23  -  26 
The  beams  shov;n  in  these  figures  have  the  parair-eters 
V   =  2200  in3 

E   =  30x10^  lb/in 

S^  =  50x10-  lb/in'' 


Figure  23.  Three  Flement  Optinum  Design  for  a 
Slrrply  Supported  Pean  under,  a  Con- 
centrated Load  Applied  at  Point  1. 
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Figure  24.   Three  Element  Optimum  Peam.  Designs  for 
Simply  SuDported  Peams  under  Concen- 
trated Loads  at  Points  2  and  3- 
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